
include("../src/qmcmary.jl")
using ..qmcmary

using Test
using LinearAlgebra

BLAS.set_num_threads(1)

"""
计算G
"""
function rawgf(allbmats, idx)
    Nt = length(allbmats)
    Nx2 = size(allbmats[1])[1]
    Rs = I(Nx2)#allbmats[1]
    for i in Base.OneTo(idx)
        Rs = allbmats[i]*Rs
    end
    Ls = Diagonal(ones(Nx2))
    for i in Nt:-1:(idx+1)
        Ls = Ls*allbmats[i]
    end
    return inv(Diagonal(ones(Nx2)) + Rs*Ls)
end


"""
计算G
"""
function rawgft(allbmats, idx)
    # g(t,0) = <c(t) c+(0)> = B(t,0)G(0,0)
    # g(0,t) = -<c+(t) c(0)> = −(1−G(0,0))B^−1(t,0) 
    g0 = rawgf(allbmats, 0)
    if idx > 0
        bt0 = prod(reverse(allbmats[1:idx]))
    else
        bt0 = I(size(allbmats[end])[1])
    end
    gt0 = bt0*g0
    ind = I(size(allbmats[end])[1])
    g0t = -(ind-g0)*inv(bt0)
    return g0t, gt0
end



"""
一次执行
"""
function run1()
    ##很可能是这里，加上exp(Matrix)导致的数值不稳定问题，导致工作不正常
    Nx = 72
    hk = lattice_chain(Nx, -1.0)
    Ui = 3*ones(Nx)
    the = rand(Nx)*pi
    nx = zeros(Nx)
    ny = sin.(the)#[0.5sqrt(2), 0.5sqrt(2)]#
    nz = cos.(the)
    #
    Nt = 200
    Ng = 5
    bidxs = Vector{Int}(undef, Ng)
    hscfg = Matrix{Int}(undef, Nt, Nx)
    for bi in Base.OneTo(Nt)
        cfg = 2(round.(rand(Nx)) .- 0.5)
        hscfg[bi, :] .= cfg
    end
    sp = default_splitting(Nt, hk, Ui; Z2=true)
    sslen, bmats, allbmats, ss = initialize_SS(Nt, Ng, sp, hscfg, nx, ny, nz)
    hscfg, bmats, allbmats, ss, gf, ph = dqmc_step(Nt, Ng, sp, hscfg, nx, ny, nz, sslen, bmats, allbmats, ss)
    sslen2, bmats2, allbmats2, ss2 = initialize_SS(Nt, Ng, sp, hscfg, nx, ny, nz)
    gf, ph = eq_green_scratch(ss)
    #gf2 = rawgf(allbmats2, 0)
    gf3, ph3 = eq_green_scratch(ss2)
    #println(gf - gf2)
    #@test all(isapprox.(gf, gf2, atol=1e-6))
    println(diag(gf))
    println(diag(gf3))
    @test all(isapprox.(gf, gf3, atol=1e-6))
    #ph2 = det(I(Nx*2) + prod(reverse(allbmats2)))
    #ph2 = ph2 / abs(ph2)
    #@test isapprox(ph, ph2, atol=1e-6)
    @test isapprox(ph, ph3, atol=1e-6)
end


run1()
